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ABSTRACT 


In simulating trajectory estimation problems, a 
rapid procedure is desirable for generating random 
sample -covariance matrices based on large numbers of 
observations. By using existing random-number gen- 
erators, an economical method is developed that yields 
* 

a matrix S whose elements have the same joint dis- 
tribution as the elements of the sample-covariance 
matrix S. 


t 



THE GENERATION OF A RANDOM SAMPLE -COVARIANCE MATRIX 


By Alan H. Feiveson 
Manned Spacecraft Center 


SUMMARY 


Trajectory estimation simulation problems make desirable a rapid procedure 
for generating random sample-covariance matrices based on large numbers of ob- 
servations. This paper first presents an algorithm for such a procedure and 
then shows its derivation from the Cochran-Fisher Theorem concerning quadratic 
forms. Finally, an example is given. 


INTRODUCTION 


In trajectory analysis, the "best" estimate of the state is a function of 

the covariance matrices R. associated with the observation stations. For 

1 

practical use, estimates must be substituted for the unknown exact R . In 
some cases, estimating the R^ directly from the observations may be desirable. 

The well-known "best”, or unbiased-maximum-likelihood-based (u.m.l.b.), 
estimator of a covariance matrix R_^ is given by 



where the are the observation vectors and n is the sample size. To sim- 

ulate a procedure where u.m.l.b. estimates are used, random matrices must be 
generated that have the same distribution as these estimates. 

* 

The obvious method of generating a matrix S , having the same distri- 
bution as S, is to generate the n observation vectors |X. ; i = 1,... n| . 

But if each vector X^ has p components, generating n observation vectors 

necessitates generating at least np i andom numbers. 



This paper presents 
using only p(p + l)/2 
np. 

■x* 

an alternate method of generating S which requires 

random numbers - usually a much smaller quantity than 


SYMBOLS 

A, A*, B, B*, C, R, W, S 

* 

, S matrices 

A i 

matrices in Cochran’s Theorem 

b. . 

i element of B 

* 

b . . 
ij 

th * 

ij element of B 

T 

C 

transpose of the matrix C 

I 

identity matrix 

i, 3, k 

indices of summation 

N(0,R) 

normally distributed with mean 0 and 
covariance matrix R 

N., N. . 
J iJ 

standardized normal random variates 

n 

sample size 

P 

size of covariance matrix (number of 
variables in one observation) 

Q 

J-l 

matrix equal to I - ^ ^ Q^. 

k=l 

«1 

T / • T 

matrix equal to y 

r . 
3 

row of matrix W 

T 
r . 

J 

transpose of r^ 
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T 


v . 
J 


w. . 


y J 


y j 


transpose of t 
random variable 

i j ^ element of W 

1 x (n - l) random vector in 
Cochran^ Theorem 

of a set of orthogonal 1 x (n - l) 
vectors 

transpose of y. 



X 2 (n - j) 

V . 

1 

0 


p x 1 vectors 

chi-square with n - j degrees of freedom 

rank of A. 

i 

p x 1 null vector 
is distributed as 


METHOD 


Let S = A/(n - l) be the u.m.l.b. estimator of a p x p covariance 
matrix R from an independent normally distributed sample of size n. It 
can be shown (ref. l) that 


A = 


n-1 

E 

k=l 


Z. Z 1 
k k 


where the p x 1 vectors ^z^.; k = 1, 2 ? . . . n - l| are independent and 
normally distributed with zero mean and covariance matrix R. 



Since R is a covariance matrix, it is semipositive definite. Therefore, 
a matrix C exists such that 


T 

CCT = R 


(3) 


It follows that the vector z, can be written 

k 


z. = Ct 
k k 


(4) 


where 


i) 


Let 


b = ( ho} 

k=l 


( 5 ) 


Then, 


CBC = C 


n-1 

z 

k=l 


T T 

Vk c =A 


( 6 ) 


Generation of A 


Let A be a generated matrix whose elements have the same joint distri- 


bution as those of A. To obtain 


* ** 
S = A 


yin - 1 ), 


it is necessary only to 


generate a matrix B 
Then, A is computed so that 


whose elements are distributed as the elements of B. 


* 

A 


* T 

CB C 


(7) 
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Hence, the problem is reduced to generating the random symmetric matrix B . 

■Xr 

An algorithm for generating B is given below. For a justification of this 
procedure, refer to the Analysis. 

Generation of B 


1. Generate p independent X variables v., j = l,...p, having n - j 

J 

degrees of freedom. One method of obtaining v^. is to generate a standard 

normal variate N and substitute it into the Wilson-Hilferty •• 2 
J 

(ref. 2). The approximation can be written 


X 


approximation 


v. « (n - j) 

J 


1 - 


9 (n - j) + N j \ 9 (n - j) 




2. Generate p(p - l)/2 independent standard normal variates N 4 . , i < 

J 

and j = 1, 2, . . .p. 


* 


3. Form the diagonal elements of B 


b . . , j = 1 , ...p| as follows : 

J J 


* 

13 ii = T 




b . . = v . + 
JJ J 


X B i/ (j ' ^ 


i=l 


* 


4. Form the off-diagonal elements of B as follows: 


b ’■ ' b ” = “’j-vp 


i,i ji "lj v i 

i-l 


b *ij ' b *n * B ij VT + 2 Vkj (i > b) 


k-1 


Once B has been generated, A follows from equation (?). 



ANALYSIS 


Using the notation of the Method section and noting that by joining the 
vectors t and k = 1, 2, ... n - 1 as columns, a p x (n-l) matrix W 

K 

can be formed 


W = { w 




where r. is the 1 x (n - l) row vector of W. Thus, the i j UIi ele- 

3 T 

ment of B, b. is equal to r.r. . 

* ij i j 


w) 


fce\ 


\^/ 


. .th 


By using the Schmidt orthogonalization process, a set of orthogonal 
vectors jy^., j = 1, 2, . . . pj. can be generated where 


y j = r j - r /i y / y i y i - ••• r j y j-i yj-i/ y j-i y j-i 


- r . ! - Ql - Q 2 - ... Q._ x 


= r.Q 


( 8 ) 


J-l 


where Q. = y. T y. j Y-Y-^ > Q = I - ^ ^ and I is the (n - l) x (n - l) 


l li/ii 

identity matrix. 


k=l 


The matrices Q, Q , . . . Q have the following significant properties: 

i <3 


1 . Q 1 , Q 2 , . . . Q ^ have a rank of one . 

2. Q.Q. = 0 for i ^ j. 

-L J 

3- Q, Q ? . . . Q are symmetric idempotents. 

i j —i 

4 . Q has rank n - j . 
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Proof 


1. The vector y\ clearly spans the entire range space of Q. . 


T T 

y . y . y . y . 

2 * Q-,-0^ = r— — FFTr = 0 because y ± y T = 0 for i ^ j 


ld KJM 


3- Clearly Q is symmetric. To show idempotence, 


T T T 

y± y^i y L n y± 

Q.Q. = . ± — £ i_ — _ = Q 

i i I r T T u 

y i y i y i y i y i y i 


and 


QQ = 


1 - Q i - Vi 


1 " °1 - ••• °i-l 


= 1-2 Q a + ... Q.^ + 


= I - | + . . . Q.j +1 | = Q 


Q i + ••• Qj-l 


4. This follows from elementary theorems on idempotent matrices (ref. 
Consider the following form of the Cochran-Fi sher Theorem. 


Theorem 


xx 


If x is a 1 x ( n - 1 ) random vector distributed N ( 0 , I ) , and i f 
k 

= ^ ^ xA^x^ the ra nk of the sum of the ' s equalling the sum of the 


i=l 


ranks of the separate A.'s is a necessary and sufficient condition for xA.x 

x 2 ~~ " ] 1 

to b e distributed as central X with degrees of freedom ( where i_s 

I TT T 

? and for xA^x , xA^x 3 ' ’ ' x ^k x to be jointly independent 

(ref. 4). 
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Note that the inner product r .r . can he written 

3 0 


r j r j T ■ r o Ir / ■ r i Q + Q 1 + • • • Q j-i r / 


o-l 


= r .Qr . + 

0 0 


Z r A r / 


k=l 


( 9 ) 


Equation (9) 
Q, Q 1 , ... 


satisfies the condition of the Theorem where the matrices 
play the role of the It therefore follows that 


T T T 

r .Qr . = r .QQ r . = r .Q 

0 0 0 0 0 


r Q 
J 


T 2 

= y i y i ~x (n - j) 

J J 


Since the y. are mutually orthogonal and normally distributed, the quantities 

y y T ( j = 1, 2, . . . p) , are mutually independent. They can be generated 
0 0 _ 2 ... . 
independently using random variables v., having the X distribution with 

n - j degrees of freedom. 


Once the set 




is given, the quantities 


a . . 
ij 


r .Q.r 
J 1 J 


1/2 


T T 

r o y i y i r o 

T 

y i y i 


1/2 




y i y i 


1/2 


( 10 ) 


being normalized linear combinations of N(0,l) variates, are themselves, 
N(0,l) variates. 

Since all the elements of the matrix ¥ are mutually independent, 

is independent of cr.,. t for 0 / 0 ’ , i < 0 3 i ’ < 0 ' • Furthermore, as a 

^■0 . 

consequence of the Theorem, it is known that for i f i', o.. is independent 
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of a.,.. Therefore, the p(p + l)/2 quantities, y.y. and a. .(j = 1, p; 

1 J J J "X J 

2 

i < j), can be generated independently, using the X random variable v. for 

J 

T 

y.y. and standardized normal variates N. for a. . • 

J J ij iJ 


■x* 

The diagonal elements of B are easily computed from equation ( 9 )* Let 


b *n - v i 


b . . = v. + 
JJ 3 


2 (d ’ 11 
i=l 


V T T 

y.y. = r.y. , it follows that 
xi j 1 


N 


ij nTT J“1 




From equation ( 7 ) tor i < j, 


r.y. = r . 
J x j 


r .r . 
J 1 


N n 


r /2 


N 2i 




r j y 2 



T) 

f T 


t\ 


T 

T* — 

( r i y i i 

T r i y 2 

T 

bi y i-l 1 

T 

y i-l 

X . — 

1 

r 

y i y i T ! 

^1 / Tl 

1 ( y 2 y 2 

y 2 " ' ' 

Tj 

7i-i y i-i 


N. , . 

X-ll 


+ ... 




r Fi-i 


ji 


N.. . N, . + N 0 .N„ . + ... N. , .N. , . 
li lj 2i 2j 1 -I 1 i-Ij 
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can be generated by 


■¥r - X - 

Therefore, b . . = b .. 

1 3 Ji 


* 


b . . = N. . . Hr — 
ij ij \ 1 


i-1 


b*. . = N. . + > N .N. .(i - 1). 

ij ij \ i Z— i ki kj' 


k=l 


Example 


Consider the generation of S based on 101 observations 



i 

ir\ 

-.21 

0 

when R ls given to be 

-.21 

.50 

• 05 


0 

-05 

.25 



r .6 

-.3 o 


Then n = 101 , P = 3 ? and C = 


0 .7 .1 

0 0.5 


It is necessary to generate only 6 (instead of 606 ) random numbers from an 
N( 0 ,l) population. They are: 


1 

^2 


= -0.258 

II 

C\J 

-0.585 

= -0.882 

OJ 

II 

0.332 

= 1.869 

N 23 = 

-0.110 
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2 

The Wilson-Hilferty X approximation gives: 


= 100 

! . 2 _ 

, (-0.238) ypT 

(9) (100) 

yf^T 

= 99 

! 2 

, (-0.882) V2 - 

( 9 ) ( 99 ) 

V89T 

= 98 

l 2 

, (-1.869) Vi" 

(9) (98) 

V882 


96.027 


86 . U92 


125.769 


Finally, the procedure given in the Method section yields 




b = 96.027 


b *22 = 86 *^92 + (-0.585) 2 = 86.835 

b* = 125.769 + (0.332) 2 + (-0.110) 2 = 125.891 


* 12 = -0.585 V96.027 = -5-73^ 


b 15 = 0.332 ^/96.027" = 3.250 


b* 2 j = -0 . 110 \f86 . 492 + (-O.585) (0.532) = -1.216 


Thus, 


* 

A 


T * 

C B C 


44.449 

-20.412 

1.157 


-20.412 1.157 

43 . 638 5 - 869 

5.869 31.473 


ll 



and 


0.444 

-0.204 

0.012 

-0.204 

0.436 

0.059 

0.012 

0.059 

0.315 


CONCLUDING REMARKS 


This report has presented an economical method of generating a p x p 
sample covariance matrix based on n observations. The method requires the 
generation of only p(p + l)/2 random numbers instead of the usually much 
larger quantity np. The matrix C referred to in the Method section may be 
obtained by methods readily adaptable to computers. 


Manned Spacecraft Center 

National Aeronautics and Space Administration 
Houston, Texas, October l8, 1965 
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